
ntl.zones <- sf::read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01') %>%
  st_set_crs(.,4326)

countries_GID_0 <- unique(ntl.zones$GID_0)

sf::st_crs(ntl.zones) <- 4326
# 102022
ntl.zones.prj <- ntl.zones %>%
  st_transform(.,crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)
ntl.zones.pts2 = dplyr::bind_cols(ntl.zones.pts,st_drop_geometry(ntl.zones))
wblkc <- read_sf(paste0(wkdir,'/local/gis_data/003_boundaries/lake_chad_latest/lake_chad_region_dissolved.shp'))

ntl.zones.pts.out <- st_join(ntl.zones.pts2[c("OBJECTID","GID_0")], wblkc, join = st_within)
out.df <- as.data.frame(st_drop_geometry(ntl.zones.pts.out))
save.dta13(out.df,paste0(wkdir,'/proc_data/WBLKC.dta'))
